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Abstract. The virial expansion of ideal quantum gases reveals some interesting and 
amusing properties when considered as a function of dimensionality d. In particular, the 
convergence radius p c {d) of the expansion is particulary large at exactly d = 3 dimensions, 
p c ( 3) = 7.1068 ... x limd-j .3 p c (d). The same phenomenon occurs in a few other special 
(non-integer) dimensions. We explain the origin of these facts, and discuss more generally 
the structure of singularities governing the asymptotic behavior of the ideal gas virial 
expansion. 


1. Introduction 

To cite an authoritative source, the treatment of Bose-Einstein and Fermi-Dirac perfect 
gases can be made in an extremely simple manner |l|. The topic has been known since the 
discovery of quantum statistics |2-5]. One might think that from a theoretical perspective 
there is absolutely nothing new to discover about it. Nonetheless, in this paper we report 
some rather amusing behavior found when investigating how the virial expansion and other 
properties depend on dimensionality d of the system. Here, d is an effective dimension given 
by d = 2D/u, where D is the dimension of physical space and v is the power of wave number 
entering the dispersion relation of the excitations in question (see below). For nonrelativistic 
massive particles, we have d = D. For massless fermions (a reasonable description of 
neutrinos and electrons in topological insulators), we have d = 6, and for the asymptotic 
low-energy spectrum of electrons in graphene, we have d = 4. Moreover Bose-Einstein 
condensates with effective spatial dimensions D ranging from 0 to 3 are now routinely 
made. Such systems are of considerable current interest. While interactions obviously play 
a role to varying degrees in the above mentioned systems, this has nevertheless motivated 
us to revisit the behavior of ideal quantum gases as a function of dimensionality with an 
emphasis on the analytical structure of the equations of state and other thermodynamical 
quantities. The analytical structure reveals itself as surprisingly intricate and complex as 
dimension is varied. 

For any dimension d the equation of state of an ideal quantum gas has a nonzero ra¬ 
dius of convergence at p = 0. It can therefore be analytically continued to a complete 
Riemann surface with much (d-dependent) mathematically interesting structure. We have 
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investigated how and where singularities occur on this Riemann surface, how they change 
character and position with changing d, and how that governs the asymptotic behavior of 
virial coefficients. We have even found cases where the analytically continued equation of 
state reappears in a new form, meaningful from the requirements that both density p and 
pressure p are real for physical values of the chemical potential p (but usually pathological 
with regard to more detailed physical behavior). To our knowledge this has not previously 
been reported in the literature. |6j 

An ideal Bose gas undergoes a condensation in which the zero-momentum ground state 
becomes macroscopically occupied above a critical density (or theoretically equivalent be¬ 
low a critical temperature) provided the dimensionality d > 2. This is a phenomenon 
accompanied by true non-analyticities in thermodynamic potentials (however without the 
appearance of collective modes). One would then expect on general grounds that any ex¬ 
pansion of physical quantities (e.g. pressure) in powers of density would have a convergence 
radius given by the critical density or smaller. 

On the other hand, it has been known for a long time that the virial expansion of an ideal 
non-relativistic Bose gas in d = 3 dimensions possesses the peculiar and somewhat surprising 
property of having a convergence radius p c { 3) which is much larger than the critical density 
for Bose-Einstein condensation, pbe(3) = C(f)Ay 3 ( w h ere A t is the thermal de Broglie 
wavelength, and ( is the Riemann zeta function). 

This was first demonstrated by Fuchs |Yj, following a conjecture by Widom 18]. Jensen 
and Henuner |9| (see also Ziff and Kincaid 0) estimated that p c ( 3) « 7 x /?be(3), based 
on numerical evaluation of the first 116 virial coefficients. The computation is by no means 
a trivial one — the n’th virial coefficient has a magnitude of order 10~ 127n , the remainder 
of cancellations between terms of order 1. A reliable evaluation of the 116’th coefficient 
required computations to be carried out to 160 digits accuracy — a quite formidable task 
in 1971. In the above works, the reason for the existence of such a surprisingly large 
convergence radius was not discussed. 

In d = 2 dimensions the situation is the opposite, pbe(2) is infinite while p c ( 2) = 27rAy 2 
(as can be seen from the repeatedly discovered exact virial expansion in two dimensions 
|10ffl3| ). This led us to investigate the dependence on dimensionality in more detail. This 
had previously been done by Ziff and Kincaid [jl0| (they focused on integer dimensions only, 
but some of their results are valid for arbitrary dimensions). Actually, our parameter d 
need not correspond to the physical dimension of the system — what matters is that the 
density of states per energy interval scales like 

sOO ~s d / 2 ~\ (1.1) 

and that it is possible to define pressure and density as position independent quantities. For 
excitations with dispersion relation e{p) ~ | p\ v in D physical dimensions we find d = 2 D jv. 
Only for a standard nonrelativistic spectrum will d correspond to the physical dimension. 
With excitations living on a fractal subset of physical space we may also obtain a non-integer 
d. Although there may be many interesting physical realizations, the purpose of varying d 
continuously in this paper is mainly to obtain a more coherent and connected picture of the 
behavior of the virial expansion as dimensionality is varied. 

The rest of this paper is organized as follows. In Section [2] we first discuss a peculiar 
symmetry between fermions and bosons which occur in the case of ideal quantum gases, its 
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origin, and how it extends to interacting systems. Next we discuss some general aspects 
of the virial expansion, and display the first few term of the ideal Bose gas expansion for 
general dimension d. 

Section [3] describes our numerical exploration of the virial coefficients A n as function 
of dimension d. Each A n has a number of zeros. They reveal intriguing patterns which 
inspired many conjectures (and eventually this whole research). 

In Section [4] we consider the special cases of d = 0 (quantum dot), d = 2 (confined layer), 
and d = 3. Apart from their physical applications, d = 0 and 2 are interesting since in 
those cases the relation between fugacity z and density p can be inverted exactly. Hence, 
they allow for a more explicit and detailed analysis, and provide boundary conditions which 
must be obeyed by the general expansions. The physically most important system, d = 3, is 
peculiar in that the equation of state (i.e., pressure p as function of p) extends analytically 
beyond the density pbe for Bose-Einstein condensation. This analytic behavior also holds 
for other thermodynamic quantities like the chemical potential (the density fluctuations 
exhibit a pole singularity at pbe)- The behavior for p > pbe is easily seen to be unphysical, 
and is not related to the correct behavior of the condensed state. 

Section [ 5 ] is a mathematical discussion of how the parametric representation (2.2) can 
be extended to the Riemann surface of the poly logarithmic functions involved, leading to 
the general parametric representation (5.7). The sheets of the complete Riemann surface 
are labeled by an infinite-dimensional vector k of integers. On an infinite subset of these 
sheets the parametric representation (5.7) provide a candidate equation of state, related to 
the usual (low-density) Bose equation of state by analytic continuation. 

Section [6] is a first account of our travel on the Riemann surface of the equation of state, 
first locating the singularities governing the radius of convergence of the virial expansion 
for d = 3, and next exploring how these singularities move as d is changed. 

In Section [7] we demonstrate how knowledge of the closest singularities of p(p), and the 
behavior of p(p) around these, can be used to provide an accurate analytical representation 
of the asymptotic behavior of the virial coefficients. The analytic prediction compares very 
well with numerically calculated coefficients. 

In Section[8]the behavior near the Bose-Einstein condensation point is analysed, revealing 
why the singularity in the equation of state vanishes at this point for dimensions d = 
2 + 2/(m + 1) (with m = 1... 6). 

Section [9] is a second account of our exploration of the Riemann surface, with focus on 
how singularities of the equation of state flow when d is changed, and in particular the 
behavior of this flow as d —> 2 and d —> 0. The singularities are in general of square root 
type. However, for d = 2 the singularities are logarithmic; hence each of them must be 
formed by an infinite number of coalescing square root singularities as d —> 2. Moreover, 
(infinitely) many pairs of square root singularities flow towards p = 0 as d —> 0, where 
none can be seen in the explicit equation of state. The annhilation mechanism seems to be 
similar to the one which is operative in the disappearance of two square root singularities 
in \/z 2 — e 2 when e —> 0. 

We close the paper with a few remarks in Section 10 
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2. VlRIAL EXPANSIONS FOR IDEAL QUANTUM GASES 

Consider a d-dimensional volume V = L d filled with identical nonrelativistic, noninter¬ 
acting mass-m quantum mechanical particles. Without internal degrees of freedom, the 
number of single-particle states in a tiny interval de around e is 


V (2nm£\ d/2 de ^ ^ ^ 

Tmv*-) T= V9{£)d£ 


( 2 . 1 ) 


when e > 0 (otherwise zero). Here, the surface of a unit sphere in d dimensions is set to 
Sd -i = 2 7r rf//2 /r(d/2) in general, and we have assumed that V — > oo in a regular manner. 
At a given temperature T we choose to measure energy in units of /c#T, and length in units 
of the thermal wavelength A t = (2'rrmkBT/h 2 )^ 1 ^ 2 . For bosons in the grand canonical 
ensemble the equation of state is given in implicit form by 

1 r°°d£Z£ 1+d / 2 

vrrnh (2 2) 

1 

P = 


e e £ — z 
r°° de z£ d ! 2 


e e £ 


r(d/2) Jo 

where z = is the fugacity and p the chemical potential. The expressions for fermions are 
the same with the replacements 


P 


-z , 


P 


-p. 


(2.3) 


It is intriguing that the same functions describe both bosons and fermions, only in different 
parameter ranges. 


However, the generalization of equation (2.3) to interacting systems is less direct. In a 


functional integral formalism, the fugacity relation between bosons and fermions is related 
to the rule that for bosons one should integrate over classical fields ip which are periodic 
in the imaginary time (t) direction, and for fermions over Grassmann fields x which are 
antiperiodic. This difference can be eliminated by a T-dependent gauge transformation at 
the cost of transforming fi —»• p + (2 n + l)vri, equivalent to z —> —z. Without interactions 
the relation between pressures p is a consequence of the relation between classical and 
Grassmann gaussian integrals, 


eP v = 


# v = 


j Vp*Vip e 1 A, 

j T>x*'E>x e _x Ax oc det A, 


which in the interacting case generalizes to the diagrammatic rule of one extra minus-sign 
per fermion loop. I.e., if we for a model with interactions know the complete loop expansion 
of the bosonic functional integral, with Pl(z) being the L-loop contribution, then relation 
(2.3) generalized to 


d 


Pf( z ) = z), Pf(z) = z—p f {z) 


where / = 0 for bosons and / = 1 for fermions. 


(2.4) 
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In the following, we return to the noninteracting model and we will mainly consider 
bosons in explicit expressions. The integrals (2.2) can be written as power series in z 
(Mayer expansions), 


p= z+ J2-^d^ = z+ J2 b ^ z ^ 


e =2 

oo 


1=2 


(2.5) 


P = z 


+ E 


£d /2 


= Z 


+ 'y~' j £bez i . 


1=2 1=2 

The fermion expressions are obtained by bi —> (— Y +l bp. By eliminating z and expressing 
p in terms of p we obtain the virial expansion, 

OO 

p = p + ^ 2 A n p n , (2.6) 

n =2 

where A n are the virial coefficients. The fermion cofficients are related to the bosonic ones by 
A n —> (—1 ) n+l A n . Therefore, the bosons and fermions have the same radius of convergence 
for their virial expansions, governed by singularities which are related by inversion, p —> —p. 
Due to reality conditions both cases also have reflection symmetry under Im p —> —Im p. 
Thus, singularities outside the real axis will occur in complex conjugate pairs. 

The virial coefficients can be calculated by order-by-order inversion of the Mayer se- 


the power series 


ries (2.5). They can also be found through an explicit algorithm |14j. For the latter, define 

(2.7) 


l=i 


g ( z ) = X> + i)W, 

S ge] 

tg(z) 

1 + tg(z) 


and coefficients C nrn through the generating function 

oo n 

- XI E CJ nrnZ n f 


( 2 . 8 ) 


n= 1 m =1 

An alternative way to describe C nm is as follows: Define v n = — (n + l)6 n+ i. Then 

n 


Cnm. — 


E'( 


V\ l>2 VZ ■ ■ ■ ) V? V^V^ 3 ■ ■ ■ 


where the sum is over all sets {u\ 122 123 ■ ■ ■ } of non-negative integers such that v k = m 
and Ylk = n - In any case the virial coefficients are 


A n +1 — 


n 

— Y 

n + 1 


m =1 


n + m — 1 
m 


c„ 


Likewise, the fugacity can be expressed in terms of density as 


z = p 


+ E Bn 


P n , 


n =2 


where 


Bn +1 — 


1 n 

1 y 

4 - i 


n + 1 


m=1 


n + m 
m 


C n 


(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 
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Dimensional dependence of low virial coefficients 



Figure 1. The 3 rd to 5 th virial coefficients as function of dimension d. 
To adjust all curves to the same plot we have multiplied A n by expressions 
e Q+,9d , with real n-dependent coefficients a and /?. Note that A n (d) has a 
zero which rapidly approaches d = 3 as n increases. 


The first virial coefficients are explicitly (with x = 1 + d/2), 

A 2 {d) = — 2 ~ x , 

Az{d) = 4 • 4~ x — 2 ■ 3~ x , 

A 4 (d ) = -20 • 8"* + 18 • 6 ~ x - 3 • 4~ x , (2.12) 

A 5 (d) = 112 • I5~ x - 144 • 12~ x + 

18 • 9~ x + 32 • 8 ~ x - 4 • 5~ x . 

They are plotted in figure [lj Fermion coefficients are obtained by A n —> (—1 ) n+ 1 A n . 


3. Numerical experiments 


The virial coefficients (2.12) exhibit much structure when analysed as function of di¬ 
mensionality. They reveal an interesting pattern of zeros. The most prominent feature is 
that one zero d\{n) rapidly approaches d = 3 as n increases, and another ^(n) rapidly 
approaches d = |. Numerically, we find 


3 — d\(n) = 2 • 10~\ 4 • 10“ 2 , 6 • 10~ 3 , 


7-10 


-4 


(3.1) 
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Figure 2. Scatter plot of the zero’s of the virial coefficients A n {d ) in the 
region 2 < d for 3 < n < 30. There are additional zeros in the region d < 2. 


Zero’s of virial coefficients A n (d) 



for n = 3,4,.... The convergence turns out to be exponential in n, with some oscillations in 
the prefactor. Clearly, it is the presence of the zeros close to d = 3 which causes the virial 
expansion to have an unusually large radius of convergence at d = 3. 

More generally, the set of zeros form an intriguing pattern at low n. Investigating this 
order by order is an amusing numerical experiment, providing many opportunities to make 
(wrong) conjectures. When extending this experiment up to n = 300, we found the most 
prominent features to be that i) at d = 2 all even virial coefficients beyond the second vanish, 
^2^(2) = 0 for n = 2,3,..., and ii) there are several special dimensions d m which act as 
“attractors” for zeros, with the consequence that the virial expansion has an exceptionally 
large convergence radius, p c (d m ) > PBc(dm)-, for such dimensions. 

We soon discovered that the special dimensions followed the pattern 

d m = 2-\ -—, m = l,2,..., (3.2) 

m + 1 


leading to the conjecture that this is true for all m. The correct result, which is difficult 
to discover by numerical experiments alone, is that equation (3.2) holds only for rn = 
1, 2,..., 6. The first 5 members of this set are clearly visible as the emerging lines in Fig. |3j 
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Figure 3. Scatter plot of the zero’s of the virial coefficients A n (d) in the 
region 2 < d, 3 < n < 300. There are additional zeros in the region d <2. 


4. Special dimensions 


We next consider some special cases of d. 
pletely when d = 0 and d = 2. 


The expressions (2.5) can be analysed com- 


4.1. Quantum dot. For d = 0 we find p = z/(l — z) and p = — ln(l — z). I.e. 

P = In (1 + p) = P + E ~^—P n i ( 4 -l) 

Z — J n 

n =2 

with virial coefficients A n = (—l) n+1 /n and radius of convergence p c { 0) = 1 due to the 
logarithmic singularity at p = —1. This singularity has no natural physical interpretation 
for bosons. However, it corresponds to the maximum obtainable density for fermions. The 
d = 0 fermionic equation of state is the same as for a classical hard core lattice gas (without 
other interactions). 

4.2. Two-dimensional layer. For d = 2, we find p = — ln(l — z) and can infer the diffe¬ 
rential equation (dp/dp) = p/(e p — 1) with p{ 0) = 0. The solution is 

= dilog(e -p ). (4.2) 

Note that dilog(e p ) = ~\p 2 ~ dilog(e _p ), as can be verified by direct manipulation of the 
integral. Thus, the equation of state for fermions and bosons only differ in the second 


P = 


dtt 
e* — 1 
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Equation of state for d — 0 



Figure 4. The equation of state for ideal quantum gases in zero dimen¬ 
sions. The boson equation of state in its unphysical third quadrant (p < 0, 
p < 0) actually describes the fermion equation of state in its physical region, 
obtained by inversion about the origin as shown. This is a general feature 
of noninteracting systems in any dimension. 


virial coefficient, Pfermion = Pboson + \p 2 ■ This can also be seen from the explicit virial 
coefficients, A n = B n -\/n\ where B n are the Bernoulli numbers. Since B n = 0 for all odd 
n > 1 the virial coefficients A n = 0 for all even n > 2 (recall that all odd coefficients are 
equal for bosons and fermions). The property that only the second virial coefficient depends 
on statistics seems to generalize to Haldane exclusion statistics |15| (interpolating between 
bosons and fermions) in two dimensions (l6) . The requirements, beyond interpretation of 
exclusion statistics, are that i) the density of states is constant as function of energy, and 
ii) a position-independent pressure and density can be defined. The latter condition is not 
fulfilled by the Hamiltonian model considered in Ref. 16, see Ref. 17, but this does not 
affect their computation of virial coefficients. 

Due to pole singularities in the integrand at t = 27rin (with integer n / 0) the function 
p(p) has logarithmic singularities at the same points. Hence, the virial expansion has a 
radius of convergence, p c { 2) = 27t. We know of no physical reasons for the occurence of 
these singularities (contrary to the d = 0 case, where the singularity has a physical origin 
in the fermion system). 

For later analysis, recall that a logarithmic singularity has infinitely many Riemann 
sheets. In this case, when p encircles the singularity at 2nim once in the clockwise direction 
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Figure 5. The equation of state for ideal quantum gases in two dimensions. 
The boson and fermion equations of state only differ in their second virial 
coefficient. Two space dimensions allow for intermediate (anyon) statistics 
which interpolates between the bosons and fermions. The equation of state 
for such hypothetical particles is not exactly known, and probably very far 
from a smooth interpolation at high densities (equivalent to low tempera¬ 
tures) . 


the function changes by p —> p — 47r 2 m. The result of repeated encirclings is that 


p —» p — 4-7T 2 ^2 m k m 


(4.3) 


when the singularity at p = 2vri m is encircled k m times in total, counted in the clock¬ 
wise direction. The order in which the encirclings occur does not matter. Hence, in the 
two-dimensional (d = 2) case, each Riemann sheet is uniquely labeled by a single inte¬ 
ger N = '^2 rn mk m . The total surface is multiply connected, so there are infinitely many 
topologically inequivalent ways of moving from one sheet to another. Amusingly, the same 


Mayer expansion (2.5) (for d = 2) defines an infinity of possible equations of state, depend¬ 
ing on which Riemann sheet is selected. All these equations of state are physical in the 
sense that the pressure is real when the density is real. 

As will be seen in Section [5] this behavior extends to all dimensions d, with the general¬ 
ization that each Riemann sheet is labeled by an infinite-dimensional vector k of integers, 
which generically must be restricted by the condition k rn = —k- m for the pressure to be 
real when the density is real. 
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Equation of state for d — 3 



Figure 6 . The equation of state for ideal quantum gases in three dimen¬ 
sions. The pressure as a function of density can be analytically continued 
beyond the critical density for Bose-Einstein condensation, pbe = C(§) = 
2.612.... The same is true for other thermodynamic quantities, like the 
chemical potential p = p(p). However, the behavior for p > /jbe is unphys¬ 
ical (with (dp/dp) < 0 and (dp/dp) < 0), and does not correspond to the 
correct infinite volume limit of the system. 


4.3. Three dimensions. In three dimensions, it seems impossible to invert the Mayer 
expansion (2.5) explicitly. We have evaluated A n (3) numerically to high accuracy and quite 
large n. They behave like 


A n (3) ~ A(n) exp (—an) cos (bn + c) 


(4.4) 


when n becomes large. Here A(n) changes quite slowly and smoothly (i.e. algebraically) with 
n, and a ~ 2.9, b ~ 0.7. As shown in Fig. [6] an equation of state can be computed from this 
virial expansion far beyond the Bose-Einstein condensation point pbe = C(§) ~ 2.612.... 
However, the solid curve is unphysical for p > pbe, since pressure decreases with increasing 
density. 

In general, the virial coefficients can be expressed by a contour integral 


A 


n 


dp_ p(p) 
2 iri p n+1 


(4.5) 


For evaluation at large n one should deform the closed curve C as far away from the origin 
as possible, since the behavior of A n then is governed by the contributions from the nearest 
singularities of p(p), beyond which C cannot be deformed. With a pair of complex conjugate 
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singularities, at p+ = p c e lu and p- = p c e 1W , one obtains an asymptotic behavior 

A n ~ p~ n cos (nw + 4>) , (4.6) 

up to algebraic corrections in n. The phase 4> is not important for locating the singularities. 
By comparing this result with the numerical observations, one concludes that the conver¬ 
gence of the d = 3 virial expansion is governed by a complex conjugate pair of singularities 
at p± ~ exp(2.9 ± 0.7i) ~ 18.5 exp(±0.7i) ~ 14± 12i. This is in agreement with earlier find¬ 
ings that the radius of convergence is much larger than the critical density for Bose-Einstein 
condensation. We want to study these singularities in more detail. 


5. Analytic continuation and Riemann surface 


There are two possible sources of singularites in an implicitly represented function p(p) = 
p(z(p)), namely i) singularities which are explicit in the parametric representation of p = 
p{z) or p = p(z), or both (in which case they may sometimes compensate each other), 
and ii) singularities occurring where p{z) is analytic, but cannot be inverted to an analytic 
function z(p) because dp(z)/dz = 0. The singularities of relevance to our case are of the 
second kind. To find their accurate positions in the p-plane, and the corresponding z-plane 


values, one must extend the parametric representation (2.2) to the full Riemann surface of 
the functions involved. 

The expressions involved are known as polylogarithmic functions, 


Li s (z) = 


1 


de ze s 


r(s) 


e e fc — z 


E 

n= 1 


(5.1) 


The sum converges for \z\ < 1 and has an analytic extension equal to the integral expression 
in the whole z-plane cut along the real z-axis from 1 to oo. This is the primary Riemann 
sheet of the polylogarithm. Moving ^ along a closed path which winds once clockwise 
around z = 1 (not encircling z = 0) changes the integral so that 

Li s (z) — t Li s (z) + — (log 2 ) s_1 , (5.2) 

T(s) 

where the right hand side again is defined on the primary Riemann sheet. The new term 
arises because moving z across the positive real e-axis drags the integration path with it. 
This can be undone by explicitly evaluating the pole contribution which makes up the 
difference. The new term must be handled carefully, since log z is multivalued around z = 0 
and (logz) 5 ” 1 is multivalued around z = 1 . We introduce p = log z as a new variable, in 
terms of which the singularity at z = 1 becomes the image of infinitely many singularities 
in the /r-plane, at 27rim, m = 0,±1,4=2,.... We define the primary Riemann sheet by 
introducing cuts 27rim + x, 0 < x < oo. If we start on the primary sheet and make k 
crossings (counted in the clockwise direction) of the cut from 27rim to oo, and no crossings 
of any other cuts, the polylogarithm will change as 


Li s (e M ) -A Li s (e M ) 

+ T(1 - s) (e 27Tik(1 - s) - l) (2vrim - p ) 3 - 1 . 


(5.3) 
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This is found most safely from the expression 18 


n=l 


^ = r o--*)(-») 


s — 1 


n =0 


(5.4) 


which is convergent in a finite region around p = 0. It can also be seen directly, in that each 
new crossing adds a new function, but also changes the phase of the previously emerged 
ones. Hence, the total contribution from k crossings is 

k -1 


27ri 


r (») „ 0 


£' 


4 — 2nin(s— 1 ) 


\s-l 


7T 


(In z 

^ e 2irik(l-s) _ ^ ( e -^i ln .) 


S— 1 


T(s) sin7r(l — s ) 

By paying close attention to phase relations on the primary sheet one finds e -7ri In z = —p. 
By further using the formula T(t)T(l—t) = n/ sin nt one finds agreement with equation (5.3). 
This equation is also valid if k is negative, that is if one makes crossings in the anticlockwise 
direction. 

The contributions from crossing different branch cuts are additive and do not interfere 
with each other. Thus, it does not matter in which order the various crossings have occurred, 
only their total number. This is a great simplification. Thus, we may label the Riemann 
sheets by an integer-valued vector k, where k m is the net number of times the branch cut 
from 27rim to oo has been crossed in the clockwise direction (starting from the primary 
sheet), and define a polylogarithnr extended to the complete Riemann surface as 


where 


Li(/r, s; k) = Li s (e M ) + ^ C(k m , s) (2nim - /r) 


C(M = j^ sinrtq-s) ^,-,) 


s— 1 


(5.5) 


(5.6) 


T(s) sin 7 r(l — s) 

When starting from the primary sheet, in practice only a few (if any) of the coefficients k m 
will be nonzero. Equation (5.3) provides the analytically continued expressions for pressure 
and density, 

p = Li(/r, 1 + d/2] k), 

(5.7) 

p = Li(^, d/2] k). 


Note that C{—k, s) = C(k , s)*. Therefore, p and p will be real when k m = —k- m for all m, 
thereby making the relation between p and p a candidate for a physical bosonic equation 
of state (with the physical region defined by p being real). Additional possibilities arise for 
rational values of s, when the logarithmic singularities at p = 27rim become algebraic and 
the coefficients C become periodic in k. Note that 

2'7ti A* 

C(k,m) = (-l) m_1 ? -— for m = 1,2,..., 

[m — lj! 

C(fc, —m) = 0 for m = 0,1, — 
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6. Locating the singularites 


It would be almost impossible to find the singularities of interest if one did not know 
where to look. Fortunately, we are in a situation quite similar to the one of explorers 
searching for a magnetic pole, in that we have a compass showing which direction to follow. 

From the numerically calculated virial coefficients for d = 3, we conclude that the behav¬ 
ior is governed by the singularity at p + « 14 + 12i and its complex conjugate. It is then easy 
to map out a path from the origin to by solving the equation p(p(t)) = tp + numerically 
and step by step for t = At, 2 At ,.... This leads to the Riemann sheet defined by Aq = 1, 
all other k m = 0. As t approaches 1, the inversion becomes difficult in accordance with our 
expectation that p(z) cannot be inverted at the singularity, namely 

^ = -Li(/i,-l + d/2;fc)=0. (6.1) 

d z z 

One may then change the search algorithm to a direct solution of this equation, which is 
explicitly 



This equation can be solved numerically to essentially aribitrary accuracy. We find 


p + 

z+ 

P+ 


= -0.322 995155 543 097 - 6.618 613 424 959 805 i, 
= 0.683 629 720405 578 - 0.238 314 132 061 880 i, 

= 14.074421676 564 36 + 12.107496 215 70789 i 
= 18.565 581330 600 061 e 0 ' 710 413678 806 621 i. 


(6.3) 


The results of Ref. 9 and 10 are in agreement with (6.3) to within the estimated accuracy. 

A z-plane map of the journey described above, to the singularity at p+, is shown in Fig. [7] 
and described in the figure caption. The other journey follows the complex conjugate path. 
As shown, one has to cross two branch cuts in the complex z-plane. When the singularity 
first has been located it is very easy to follow its movement as one changes the dimension d. 
We discuss this in more detail in the captions to Fig. [7] and Fig. [8j and further in Section |9| 


7. Asymptotic virial coefficients 

With accurate information on the singular behavior we may use ( |4.5[ ) to compute the 
asymptotic behavior of A n as n —> oo. We may expand p{p) and p{p) around p + , that is, 
with A/i = p, — p + , Ap = p + — p, and A p = p — p + , where 

A ' 0 = “E r 7 a ^’ a p = EIf(7.1) 

' n\ n\ 

n=2 n=l 

Here, p n = r n _i = Li (p+, 1 — n + d/2; k). In particular, p\ = p + and p2 = r\ = 0 (the 
singularity condition). We may now express Ap order by order in A p, thus 

A/i = (-2/r 2 ) 1/2 y/A~p H- , 

Ap = (-2/r 2 ) 1/2 p+ \J~Ap ^ - 


(7.2) 
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Figure 7. This figure maps our explorations of singularities. The fulldrawn 
line (colored red) is the first branch cut, running from 2 = 1 to z = 00 along 
the positive real axis. The dashed line (colored blue) is the second branch 
cut, running from z = —00 to z = 0 along the negative real axis. We first 
found the (fulldrawn black) line in the z-plane which corresponds to the 
curve p(t) = t p+(3), 0 < t < 1, starting at z = 0. The filled circles (colored 
green) mark the progression as t increases. As can be seen, one must cross 
two branch cuts before reaching the singularity p+(3). 

Having found /0+(3) and equation (6.1) determining p+(d), it is easy to 
explore how the singularity moves when one changes d. This is shown by the 
next fulldrawn (black) line. The black dots and filled circles (colored red) 
mark the progression as d decreases. When d —> 0 the singularity moves to 
z = 00 (multiplied by a phase), corresponding to a singular point at p = — 1. 
Actually, (6.1) becomes ambiguous when d —> 2, 0, —2,.... To handle this, 
we may write d = 2 + ee 1¥3 with £ a small positive number (0.001), and varied 
<p from 0 to 7T. In the z-plane the singularity moves slightly above z = 0 and 
back across the second branch cut. Alternatively, one may vary p from 0 to 
an infinite number of other possible endpoints (2 k + 1)tt. Each choice leads 
to a different solution. There are infinitely many singularites approching 
p = 27ri as d approaches 2. We have mapped out several endpoint choices 
before decreasing d to 0 or increasing d back to 3. The two black dotted 
lines show how the singularity moves when the encircling is chosen to ±27 t, 
and d is increased back to 3. See the caption to Fig. [8] for a further descrip¬ 
tion. 
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Figure 8. This figure shows how the singularities p+(d ) and /Obe(gQ move 
as one changes the dimension d. Starting near d = 3 the singularity Pbe (d) 
is closer to the origin than p+(d), thus determining the convergence radius of 
the virial expansion (except when d is equal to one of the special points d m ). 
However, as d decreases p+{d ) moves inwards while pbe(^) moves outwards. 
Thus, there is a crossover dimension d c , indicated by the filled circles (colored 
green), where they are equally far from the origin (as indicated by the green 
dashed circle). We find d c « 2.252 563 996. For d < d c the convergence 
radius of the virial expansion is always determined by p+(d) (and its complex 
conjugate p-(d)). Since 18/8 < d c < 16/7 the sequence of special points d s 
ends at d s = 16/7. As we decrease d further towards 2 one discovers that 
d = 2 is a logarithmic singular point for p+(d). Hence, one may encircle 
d = 2 in various ways before changing d further. The figure displays the two 
behaviors when d is decreased to 0 (after encircling d = 2 by angles ±7r). In 
one case the singularity moves to p = — 1 when d —> 0, which is a logarithmic 
singularity of the d = 0 equation of state. In the other case p + moves to 
p = 0 when d —> 0. There is no singularity at p = 0 in the d = 0 equation of 
state, but it is possible for two square root singularities to annihilate when 
they coalesce. Finally, the figure shows the behavior when d = 2 is encircled 
by 2-7T, and increased back to d = 3. In this case p+ move to p = 0, where it 
appears to annihilate with p _. Thus, the analytically continued equation of 
state has singularities close to p = 0 when d is close to 3 and 0, but not for 
d exactly 3 or 0. 
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m 

0>m 

OLrn 

0 

0.489 092 994 674 599 

-0.115 262 558 466 782 

1 

0.807395 011652 565 

-0.639 351821397655 

2 

3.094 390113185 892 

-1.159 038 217120 639 

3 

17.394 603 934 793184 

-1.466 809 677404 770 


Table 1. Numerical parameters for the asymptotic 
virial coefficients A n in 3 dimensions. Equation (7.3), 


expansion of 
with |/cq_| = 


exp (2.921 309 400 394924) and ui = 0.710 413 678 806 621, provide an accu¬ 
rate repr esen tation for large n. The relative error in A n when summing the 
series in (7.3) to m = k can be expected to be of order n ^( fc + 1 )(a fc+1 /ao)• 


We insert this expansion for the pressure into the integral expression (4.5), and deform 
the integration contour around the square root branch cut starting at p = p+. Writing 
p = p+(l + t) gives a contribution 

■e'—PAAW f-" rj^+... 

Jo 


7r 


(1 + t) n+1 


\/2 t r 2 


3 / 2 —n _ 3/2 , 

p+ n ' H- 


to (4.5). There is a similar contribution An ' = An from the singularity at p_. There 


are also higher order contributions from weaker singularities (proportional to Ap m+1 / 2 ) in 
A p. This leads to the following expansions 


A n S3 |p+| 3/2 71 ^2 Cm COS (u> 71 + <j> m ) / 
m— n •'O 


00 dtt 1/2+rl 


= n -3/2 U. |3/2 — 


(1-M) n+1 
\p+\ 3/2 ~ n ^2 n ~ m a m cos (w n + a m ), 


(7.3) 


with amplitudes a m , c m and phases a m , <t> m which are straightforward to compute. The 
first numerical values of a m and a m for dimension d = 3 are listed in table [1} 

Because of the oscillatory behavior, it is difficult to find an accurate asymptotic fit to 
A n directly from the numerical series, in particular if the prefactor ?r^ 3 / 2 is unknown. 
In retrospect one realizes that this n-dependence is the generic behavior of the algebraic 


prefactor, due to the square root type of the singularity (7.2). 


8. Behavior at the Bose-Einstein condensation point 


To understand why the dimensions d m are special for the virial expansion, with a radius 
of convergence much larger than peE; we must investigate the equation of state for p near 
Pbe- The formula (5.4) is useful for this. We restrict analysis to the interval 2 < d < 4. It 
follows that 


(pBE ~ p) = -r(l - ^d) (-/i) 1+d/2 + ■ • • , 
(pbe-p) = C(t|) (~P) + ■■■ , 
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Large n expansion of scaled virial coefficients 



Figure 9. The exactly computed values of scaled virial coefficients com¬ 


pared with the expansion (7.3). The expansion does not work well for low 


n, due to contributions from additional singular points (further away from 
the origin) of the equation of state. 


as p —> 0 
state as p 


with pbe = C(|) an d Pbe = C(1 + %)■ 
» Pbe from below: 


Elimination of p gives the equation of 


(pbe - p) = Pbe 


(Pbe ~ p) 

-n-m 


+ 


( 8 . 1 ) 


where 5 = 2/(d — 2). This term is non-singular when 5 is equal to an integer 1 + m > 1. 
This corresponds to dimensions d = 2 + 2/(m + 1) with m = 1,2,..., giving pbe — p ~ 
(pbe — p) m+1 ■ For odd m this equation of state is obviously unphysical when p > pb, 
but this is not equally obvious when m is even. There is, however, a signature of singular 
behavior in the density fluctuations, which diverge like (pbe — p)~ m ■ 

One must al so ch eck that there are no singularities to higher orders. For this it follows 
from equation (5.4) that p and p are analytic functions of £ = (— p) 1 ^ in a region around 
p = 0 (the power series converges), with dp/d£ / 0 at £ = 0. Therefore, the function p(£) 
can be inverted. This implies that £, p = — and p are analytic functions of p in some 


region around pbe for all special dimensions given by equation (3.2). However, this does 


not necessarily imply that there will be a sudden increase of the convergence radius of the 
virial expansion at these dimensions, since the convergence may be governed by singularities 
closer to the origin than pbe- 
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R(n,D) 



1/n 


Figure 10. The ratios R(n, d) = A n (d) / A EE (d) as function of n -1 , plotted 
for n = 50, ...,500. The A n (dy s are calculated numerically (using 750 
decimals accuracy), and divided by the leading order contribution (8.2) to 
A EE (d). The first correction to this ratio is of order nr 1 when 5 > 1 (i.e. 
for d < 3), and of order n 1-5 when 5 < 1. It is straightforward to calculate 
such corrections, but they are best read by computers. 


Returning to general dimensions, 2 < d < 4, we insert equation (8.1) into (4.5) to obtain 


the asymptotic contribution to the virial coefficients from the Bose-Einstein singularity. By 
deforming the integration contour around the branch cut starting at p = pbe and performing 
the integration, we find 


4 BE ) (d) = 13(d) B(l + 6,n-5) + 


( 8 . 2 ) 


as n -) oo. Here B(l + S,n — 5) is the beta function, it behaves like n ( 1+<5 1 as n—> oo. The 
coefficient in front is 


1 


13(d) = — [— T(1 — d/2)] 5 sin(7r<5). 


7T 


(8.3) 


The Gamma function has a pole singularity as d —> 4 (i.e. <5 —>■ 1), cancelling the zero in 
sin7r<5. Thus, the equation of state is nonanalytic at />be for d = 4 (which corresponds to 
m = 0 in equation (3.2)). This can also be seen by direct analysis of the equation of state 
for d = 4. 


We have checked the accuracy of (8.2) for a few non-special dimensions, see figure 10 The 
higher order corrections have the form of a double series in powers of nY 1 and n 1_ " Since 
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Figure 11. This figure illustrates how a few (out of an infinite number) of 
singularities flow away from p = 27ri as d is increased from 2 to 3. There is 
a similar flow related by complex conjugation and ko = 1 —> ko —> — 1. For 
exactly d = 3 the latter distinction vanishes: The generically logarithmic 
singularities at p, = 27rim become square root singularities at exactly d = 3. 
As a consequence a translation symmetry under k —> k + 2 emerges. For this 
reason the apparent singularity at p = 0 for d = 3 is probabably annihilated 
by its complex conjugate, and thus not present at d = 3. 


5 —> 1 + as d —> the convergence towards (8.2) becomes slow near d = 4, as demonstrated 

by the d = 3.5 case in Fig. IT 

It is now easy to understand why A n (d ) = 0 for d « d m . The contribution from (8.2) 
vanishes like (d — d m ) It can be cancelled by the contribution from the singularities at 
p±. The latter behaves like |p±| -n cos(nu; + (j>). Hence (ignoring algebraic prefactors) there 
will be a zero when 

' Pbe' 

]p± I 


d — dn 


cos(ncj + 4>). 


(8.4) 


9. Exploring dimensionality 

For general dimension p + (and similar singularities) is determined by equation ( |6.1[ ). 
Since we know a complex conjugate pair of solutions for d = 3, defined by (p + ,ko = 1) 
and (p- = p*i_,ko = —1), it is straightforward to explore how they move with changing 
dimension, by changing d in small steps and changing k when branch cuts are crossed. A 
map of this exploration (in the z-plane) is shown in Fig. [7J and described in the last part 
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Flow of ^-singularities as d decreases from 2 to 0 



Figure 12. This figure illustrates how a few (out of an infinite number) 
of singularities flow away from p = 2-7ri as d is decreased from 2 to 0. There 
is a similar flow related by complex conjugation. 


of its caption. As d = 2 the singularity approaches z = 0, which in this case corresponds to 
p = 27ri. The latter can be seen in Fig. [8| where a map of the same curve is shown in the 
p-plane. 

Actually, when approaching d = 2 one encounters a puzzle. The singularity p+(d) is 
generally of square root type, with two Riemann sheets attached locally. However, for 
d = 2 it is logarithmic with infinitely many attached Riemann sheets. How can a square 
root singularity suddenly turn into a logarithmic one? The answer is that this is impossible 
for a single square root, but it may be possible if we have infinitely many of them. Hence, 
one must conclude that there are infinitely many p-plane singularities approaching 2-/rin 
when d —> 2. They are actually easy to locate by direct analysis of equation (6.1). They 
correspond to Re p ~ log \ d — 2|, and values of Rn p which change in steps of approximately 
2tt. 

A simple way to generate them numerically is to follow the solution as d encircles d = 2 
in the plane of complex dimensions, d = 2 + with some small positive 5. Note that 


equation (6.1) has only pole singularities as function of d, hence each circle leads back to the 


same equation. However, each change of cf) by 27r leads to a new solution. Having generated 
many new singularities this way one may again follow their paths back to d = 3. This is 
shown in Fig. US 

Similarly, we may change <j) to <j) + (2 n + 1)tt for various n, and see how the singularities 
flow when we decrease d from 2 — 5 towards d = 0. This is illustrated in Fig]12| To 


our 
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initial surprise, most singularities (probably infinitely many — all but one) flow towards 
p = 0. On the other hand, we know from the explicitly found equation of state at d = 0 
that there is no singularity at p = 0 on any Riemann sheet. Hence, it must be that they 
annihilate at d = 0. By direct analysis of the equation it is possible to see that there must 
be an infinite number of solutions to (6.1) approaching p = 0 as d —> 0. The analysis is 
similar to the one for d = 2. We have not investigated the “annihilation process” in detail. 
Another surprise is that only one square root singularity flows toward the logarithmic one 
at p = —1 for d = 0 (there is one more related by complex conjugation). The infinity of 
additional singularities needed to build a logarithmic one must flow from the other d = 2 
singular points 27rim. 


10. Summary and Conclusions 


In this paper, we have investigated the analytical structure of the virial expansion and 
equation of state of ideal quantum gases in arbitrary (real-valued) dimensions, with an 
emphasis on locating the singularities that determine the radius of convergence of the ex¬ 
pansion. These simple systems have a surprisingly rich complex analytical structure. 

When investigating the behavior near d = 3, one finds that we live in a very special 
dimension from the point of view of the virial expansion. Namely, for d close to 3 the equality 
p c {d) = PBc{d) holds for every d / 3, but not for exactly d = 3. The virial coefficients 
A n (d) have zeros which approach d = 3 very rapidly as n increases. This leads to a sudden 
increase of the convergence radius at exactly 3 dimensions, paradoxically extending the 
radius of convergence far beyond the critical density for Bose-Einstein condensation (in the 
bosonic case). Enlarging the investigations to a wider d-range, one finds the same behavior 
at a few other special dimensions d m . There are six such cases in all, d m = 2 + 2/(m + 1) 
for m = 1, 2 ... 6. We have revealed a rich analytical structure in the virial expansion, when 
suitably extended to the complex plane, that explains this behavior. 

When searching for the complex conjugate pairs of singularities which determine p c (d m ), 
we find that they are determined by the simple equation (6.2) in d = 3 dimensions, and 


more generally by equation (6.1) 


While we have focused exclusively on noninteracting quantum gases in this paper, we 
expect that some of the results may be extended to interacting cases. Interacting systems 
where mean-field theories are applicable (that is, above some lower critical dimension) are 
essentially effective single particle problems and as such may fall within the class of problems 
we have considered here, provided that the resulting mean-field theory features long-lived 
excitations above the ground state condensate which are either fermionic or bosonic. Mean 
field theories go beyond any order in perturbation theory and often capture interesting 
physics of strong-coupling fixed points. They tend to be exact as the dimension approaches 
an upper critical dimension. In this context, we note here that the effective dimension d 
is twice the dimension of physical space in the relativistic cases. Low-energy excitations of 
quantum spin systems defined on fractal lattices may also be described by dilute Bose gases 
with noninteger d. 

A.S. was supported by the Research Council of Norway, through Grants 205591/V20 and 
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